NASA  Contractor  Report  198279 
ICASE  Report  No.  96-8 


ICASE 


MULTI  DIMENSIONAL  ASYMPTOTICALLY 
STABLE  4TH-ORDER  ACCURATE  SCHEMES 
FOR  THE  DIFFUSION  EQUATION 


Saul  Abarbanel 
Adi  Ditkowski 


’'Apijrowi  te 

,|>lsrtiiboiaEm  PwacitM  j 


NASA  Contract  No.  NAS 1-19480 
February  1996 

Institute  for  Computer  Applications  in  Science  and  Engineering 
NASA  Langley  Research  Center 
Hampton,  VA  23681-0001 

Operated  by  Universities  Space  Research  Association 


National  Aeronautics  and 
Space  Administration 

Langley  Research  Center 
Hampton,  Virginia  23681-0001 


®llc  quality  iuspectpju  ar 

1 996041 1  1 44 


Multi-dimensional  asymptotically  stable  4^^-order 
accurate  schemes  for  the  diffusion  equation. 

Saul  Abarbanel* 

Adi  Ditkowski* 

School  of  Mathematical  Sciences 
Department  of  Applied  Mathematics 
Tel- Aviv  University 
Tel- Aviv,  ISRAEL 


Abstract 

An  algorithm  is  presented  which  solves  the  multi- dimensional  diffusion  equation 
on  complex  shapes  to  4*^-order  accuracy  and  is  asymptotically  stable  in  time.  This 
bounded-error  result  is  achieved  by  constructing,  on  a  rectangular  grid,  a  differentiation 
matrix  whose  symmetric  part  is  negative  definite.  The  differentiation  matrix  accounts 
for  the  Dirichlet  boundary  condition  by  imposing  penalty  like  terms. 

Numerical  examples  in  2-D  show  that  the  method  is  effective  even  where  standard 
schemes,  stable  by  traditional  definitions,  fail. 
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1  Introduction 


Recently  there  has  been  renewed  interest  in  finite-difference  algorithms  of  high  order  of 
accuracy  (4*^  and  above),  both  for  hyperbolic  and  parabolic  p.d.e’s  (see  for  example,  [1],  [2], 
[3]  ).  The  advantages  of  high-order  accuracy  schemes,  especially  for  truly  time  dependent 
problems,  are  often  offset  by  the  difficulty  of  imposing  stable  boundary  conditions.  Even 
when  the  scheme  is  shown  to  be  G.K.S.-stable  the  error  may  increase  exponentially  in  time. 

This  paper  is  concerned  with  4‘**-order  approximations  to  the  long  time  solutions  of  the 
diffusion  equation  in  one  and  two  dimensions,  on  irregular  domains.  By  an  irregular  domain, 
we  mean  a  body  whose  boundary  points  do  not  coincide  with  nodes  of  a  rectangular  mesh. 

In  section  2  we  develop  the  theory  for  the  one-dimensional  semi-discrete  system  resulting 
from  the  spatial  differentiation  used  in  the  finite  difference  algorithm.  Energy  methods  are 
used  in  conjunction  with  “SAT”  type  terms  (see  [1]),  in  order  to  find  boundary  conditions 
that  preserve  the  accuracy  of  the  scheme  while  constraining  an  energy  norm  of  the  error  to 
be  temporally  bounded  for  alH  >  0  by  a  constant  proportional  to  the  truncation  error. 

In  section  3  it  is  shown  how  the  methodology  developed  in  section  2  is  used  as  a  building 
block  for  the  multi- dimensional  algorithm,  even  for  irregular  shapes  containing  “holes.” 

Section  4  presents  numerical  results  in  two  space  dimensions  illustrating  the  long-time 
temporal  stability  of  the  method,  in  contradistinction  to  “standard”  methods  for  cartesian 
grid  on  irregular  shapes. 
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2  The  One  Dimensional  Case 


We  consider  the  following  problem 


du 

dt 

=  A:— +  /(x,t);  Ti,  <  X  <  Tfi,  t  >  0,  A:  >  0 

(2.1a) 

u(x,0)  =  Uo(x) 

(2.1b) 

u(rL,t)  =  giit) 

(2.1c) 

II 

(2.1d) 

and  C^. 

Let  us  spatially  discritize  (2.1a)  on  the  following  uniform  grid: 


^2  ^5  ^N.2 

Figure  1:  One  dimensional  grid. 


Note  that  the  boundary  points  do  not  necessarily  coincide  with  Xi  and  xj^.  Set  Xj^i  —  Xj  =  h, 
1  £  i  <  —  1;  3;i  —  Tl  —  jLh,  0  <  7£,  <  1;  Tr  —  x^  =  'ynh,  0  <  'Jr  <  1. 

The  projection  unto  the  above  grid  of  the  exact  solution  u{x,t)  to  (2.1),  is  Uj{t)  = 
u{xj,t)  =  u(t).  Let  .D  be  a  matrix  representing  the  second  partial  derivative  with  respect  to 
X,  at  internal  points  without  specifying  yet  how  it  is  being  built.  Then  we  may  write 

— u(t)  =  fe[^'u(t)  +  B  +  T]  +  f(t)  (2-2) 


2 


where  T  is  the  truncation  error  due  to  the  numerical  differentiation  and  f(t)  =  f(xj,t), 
1  <  i  <  The  boundary  vector  B  has  entries  whose  values  depend  on  gL,9R,  'ILiIr  in 
such  a  way  that  +  B  represents  the  2“*^  derivative  everywhere  to  the  desired  accuracy. 
The  standard  way  of  finding  a  numerical  approximate  solution  to  (2.1)  is  to  omit  T  from 
(2.2)  and  solve 

jv{t)  =  k{bw{t)  +  B)  +  f(t)  (2.3) 

where  v(i)  is  the  numerical  approximation  to  the  projection  u(i).  An  equation  for  the 
solution  error  vector,  e(t)  =  u{t)  -  v(t),  can  be  found  by  subtracting  (2.3)  from  (2.2): 

=  kbe{t)  +  kT{t)  (2.4) 

Our  requirement  for  temporal  stability  is  that  ||  e  ||,  the  L2  norm  of  t,  be  bounded  by  a 
“constant”  proportional  to  h”^  (m  being  the  spatial  order  of  accuracy)  for  all  t  <  00.  Note 
that  this  definition  is  more  severe  than  either  the  G.K.S.  stability  criterion  [4]  or  the  definition 
in  [1]. 

It  can  be  shown  that  if  D  is  constructed  in  a  standard  manner,  i.e.,  the  numerical  second 
derivative  is  symmetric  away  from  the  boundaries,  and  near  the  boundaries  one  uses  non 
symmetric  differentiation,  then  there  are  ranges  of  values  of  7^  and  7x,  for  which  D  is 
not  negative  definite.  Since  in  the  multi-dimensional  case  one  may  encounter  all  values  of 
0  <  1L-,1R  <  I5  this  is  unacceptable. 

The  rest  of  this  section  is  devoted  to  the  construction  of  a  scheme  of  4*^  order  spatial 
accuracy,  which  is  temporally  stable  for  all  jiy'fR. 
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The  basic  idea  is  to  use  a  penalty-like  term  as  in  the  SAT  procedure  of  ref  [1];  here, 
however,  it  will  be  modified  and  applied  in  a  different  manner. 

Note  first  that  the  solution  projection  Uj{t)  satisfies,  besides  (2.2),  the  following  differ¬ 
ential  equation: 

^  =  kDn  -b  fcTe  -1-  f(t)  (2.5) 

at 

where  now  D  is  indeed  a  differentiation  matrix,  that  does  not  use  the  boundary  values,  and 
therefore  Te  7^  T  but  it  too  is  a  truncation  error  due  to  differentiation. 

Next  let  the  semi-discrete  problem  for  v(i)  be,  instead  of  (2.3), 

^  =  k[Dv  -  tl{Aly  -  Zl)  -  7‘b(Abv  -  g/?)]  +  f(t)  (2.6) 

at 

where  =  (1,  ■  •  • ,  1)^5'l(0;  Sh  =  (1,  •  •  • ,  are  vectors  created  from  the  left  and 

right  boundary  values  as  shown.  The  matrices  and  Ar  are  defined  by  the  relations: 

Alvl  =  gL-  Tl;  Ai?u  =  gR-  Tr,  (2.7) 

i.e.,  each  row  in  Al{Ar)  is  composed  of  the  coefficients  extrapolating  u  to  its  boundary  value 
gLigR.),  at  r£,(r/i)  to  within  the  desired  order  of  accuracy.  (The  error  is  then  Tl(T/j).)  The 
diagonal  matrices  tr  and  tr  are  given  by 

Tr  =  diag  {tr, , ri,2 . . . , tl^);  tr=  diag  {tr, ,...tr^)  (2.8) 

Subtracting  (2.6)  from  (2.5)  we  get 

—  k[De  —  trArZ—  trArZ d-  Ti]  (2.9) 

at 
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where 


Ti  =  Te  +  TiTi  +  TrTr 

Taking  the  scalar  product  of  e  with  (2.9)  one  gets: 

1  d 

2di  ^  ~  ~  ”^1) 

=  k{t, Mt) -\- k{e,  Til)  (2.10) 

We  notice  that  (e,  Me)  is  (e,  (M  +  M^)e72,  where 

M  =  D  -  tlAl  -  trAr.  (2.11) 

If  M  +  M'^  can  be  made  negative  definite  then 

(e,  (M  +  M^)e/2  <  -co  ||  e  |i^  (co  >  0).  (2.12) 

Equation  (2.10)  then  becomes 

5^  l|e'P<  -ic„||?|P  +i(?,T.) 

and  using  Schwartz’s  inequality  we  get  after  dividing  by  ||  e  || 

^l|e-||<-ic<,||e||+fc||T.  II 

and  therefore  (using  the  fact  that  v(0)  =  u(0)) 

||r||<ilL!!«(i_e-‘=‘»<)  (2.13) 

Co 
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where  the  ^constant"  ||  Ti  ||a/=  niaxo<T<t  ||  Ti(t)  ||. 

If  we  indeed  succeed  in  constructing  M  such  that  M  +  is  negative  definite,  with  cq  >  0 
independent  of  the  size  of  the  matrix  M  as  it  increases,  then  it  follows  from  (2.13)  that 
the  norm  of  the  error  will  be  bounded  for  all  t  by  a  constant  which  is  0{h'^)  where  m  is 
the  spatial  accuracy  of  the  finite  difference  scheme  (2.6).  The  numerical  solution  is  then 
temporally  stable. 

The  rest  of  this  section  is  devoted  to  this  task  for  the  case  of  m  =  4,  i.e,  a  fourth  order 
accurate  finite  difference  algorithm. 

Let  the  n  x  n  differentiation  matrix,  D,  be  given  by 


45 

-154 

214 

-156 

61 

-10 

10 

-15 

-4 

14 

-6 

1 

-1 

16 

-30 

16 

-1 

-1 

16 

-30 

16 

-1 

-1 

16 

-30 

16 

I2h^ 


-1 

16 

-30 

16 

-1 

-1 

16 

-30 

16 

-1 

1 

-6 

14 

-4 

-15 

10 

-10 

61 

-156 

214 

-154 

45 

(2.14) 

The  upper  two  rows  and  the  lower  two  rows  represent  non-symmetric  fourth  order  accurate 
approximation  to  the  second  derivative  without  using  boundary  values.  The  internal  rows 
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are  symmetric  and  represent  central  differencing  approximation  to  u^x  to  the  same  order. 
Note  that  D  is  not  negative  definite,  and  neither  is  the  symmetric  part  of  |(D  +  D'^)  which 
is  given  by: 


90  -144  213  -156  61  -10 

-144  -30  12  13  -6  1 

213  12  -60  32  -2  0 

-156  13  32  -60  32  -2 

61  -6  -2  32  -60  32  -2 

-10  1  0  -2  32  -60  32  -2 

1 

m2 

-2  32  -60  32  -2  0  1  -10 

-2  32  -60  32  -2  -6  61 

-2  32  -60  32  13  -156 

0  -2  32  -60  12  213 

1  -6  13  12  -30  -144 

-10  61  -156  213  -144  90 

(2.15) 

In  order  to  construct  M  we  need  to  specify  Al,  Ar,  tr  and  tr.  We  construct  A^  as 
follows: 

Al  =  +  CRAf^  (2.16) 
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where 


’  oil 

Ot2 

0(3 

0:4 

OC5 

0  . 

Oil 

0(2 

0(3 

a4 

Ol5 

0  . 

oci 

Qf2 

OC3 

a4 

OC5 

0  . 

oci 

Oi2 

Oi3 

a4 

015 

0  . 

Oil 

Oi2 

Oi3 

a4 

as 

0  . 

. 

OC2 

OC3 

Qf4 

as 

0  . 

(2.17) 


CL  =  diag  [-20ai/71,  0, . . . ,  0] 


■  -1 

5 

-10 

10 

-5 

1 

0 

...  0  ■ 

II 

-1 

5 

-10 

10 

-5 

1 

0 

...  0 

.  -1 

5 

-10 

10 

-5 

1 

0 

...  0 

The  a’s  are  given  by 


(2.18) 


(2.19) 


^2  =  -  (4u  +  f  tI  +  57!  +  571) 


IQ  1 

«3  =  371,  + -7l  +  27£  + 


(2.20) 


To  To  1  A 

-4  =  -(-7.  +  -72  +  g7!+g7lj 
=  ^-  +  ^7l  +  j7!  +  ^7l 
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Note  that  gives  a  vector  whose  components  are  the  extrapolated  value  of  v  at  a:  =  Fx, 

(i.e.,  uri(O),  to  fifth  order  accuracy;  while  gives  a  vector  whose  components  represents 

{d^v\l dx^)h^ .  Since  Cl  (see  2.18)  is  of  order  unity,  then  Alv  =  {A^^^  +  clA^^^)'v  represents 
an  extrapolation  of  v  to  to  fifth  order. 

Before  using  Al  in  (2.11)  or  (2.6)  we  must  define  tl'. 

tl  =  ^^2,  Ta,  r4,  rg,  0, . . . ,  0]  (2.21) 

where 

Tl  =  71  /2q;i 

T2  =  (-94  -  Qr2Ti)/ai 

Tz  =  (113  —  Q;3Ti)/ai  (2.22) 

r4  =  (— 56  —  Q;4Ti)/ai 

Ts  =  (11  -  OisTx)lai 


The  right  boundary  treatment  is  constructed  in  a  similar  fashion,  and  the  formulcie  corre¬ 


sponding  to  (2.16)  -  (2.22)  become: 

Ar  =  d-  CrA^^, 


■  0 . 

...  0 

0 

OCJ^-A 

OCN-2 

OCN-1 

OCN 

0 . 

...  0 

0 

OLN-i 

(XN-3 

OCN-2 

(Xn-I 

OCN 

0 . 

...  0 

0 

Oi-N-A 

^‘N-S 

OLN-2 

OCN-l 

OCN 

0 . 

...  0 

0 

OlN-4 

«iV-3 

OCN-2 

OCN-1 

OCN 

0 . 

.  ...  0 

0 

ajv-4 

O^N-3 

OiN-2 

OCN-1 

OCN 

0 . 

.  ...  0 

0 

OiN-4 

OiN-3 

OCN-2 

OCN-1 

OCN 

0 . 

.  ...  0 

0 

aN-4 

OiN-3 

OCN-2 

OCN-1 

OCN 

.  0 . 

.  ...  0 

0 

O‘N-4 

(XN-3 

OCN-2 

OCN-1 

OCN  . 

(2.23) 


(2.24) 
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(7r  =  diag[0, 0, . . . ,  0,  -20Q;iv/71] 


(2.25) 


■  0 

0 

...  0  1 

-5 

10 

-10 

5 

-1  ■ 

4«)  = 

0 

0 

...  0  1 

-5 

10 

-10 

5 

-1 

(2.26) 

0 

0 

...  0  1 

-5 

10 

-10 

5 

-1 

The  a’s  are  here: 


,  25  35  2  53  1  4 

CM  =  1  +  Y27R  +  +  jjTs 


OCN-1 


( .  13  3  3  3  1  4 

=  -  {^47fi  +  + -7j5  + -7k 


IQ  1 

«iV-2  =  37K  + J7U27I+47A 


(2.27) 


QJiV-S  = 


4  7  2  73  1  4\ 

+  gTjl  +  qIr  +  g7fi  J 


1  11  2  1  3  1  4 

aAr-4  =  j7«  +  ^7h  +  47k  +  ^7k, 


TR  —  22^2  •  *  •  i'^N-4,TN-3yTN-2,TN-l,  'Tv], 


(2.28) 


Tjsr  =  71/2Q;jv 
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tn-1  =  {—^^  —  o:n-iTn)/ocn 

TN-2  —  {IIS  —  aN-2TN)/oiN 
TN-S  =  —  OlN-zTN)/aN 


(2.29) 


TN-4  =  (11  —  «JV-4TAr)/Q^iV 
We  are  now  ready  to  construct 

i(Jl/  +  M^)  =  l{D  +  0^  -  [Ti(4‘)  +  ci4‘))  +  TK(4*>  +  c„4*>)I 

-  [ti,(4^>  +  ciAi‘‘^)  +  tb(A<,'’)  +  cjiAf »)]^}  (2.30) 

Upon  using  equations  (2.14)-(2.29)  in  (2.30)  one  gets: 


M  +  M^ 
2 


0 


0 

0 

-2  0 

32  -2  0 

-2  32  -60  32  -2 

-2  32  -60  32  -2 

-2  32  -60  32  -2 

-2  32  -60  32  -2 

0  -2  32 

-2 

0 

0 


(2.31) 
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where  and  are  6  x  6  blocks  given  by: 


Wi^)  =  (2.32) 

W(R)  =  wf )  +  )  (2,33) 


0  e  =  1  or  y  =  1 


(^OC{Tj  +  QIjTj)  tjj  ^  1 


1  <  i,j  <  5 


(2.34) 


■ -1  0  0  0  0  O' 

0  -30  12  13  -6  1 

y^(L)  ^  0  12  -60  32  -2  0 

^  0  13  32  -60  32  -2 

0  -6  -2  32  -60  32 

.0  1  0  -2  32  -60  . 

■  -60  32  -2  0  1  0' 

32  -60  32  -2  -6  0 

y^(R)  ^  -2  32  -60  32  13  0 

^  0  -2  32  -60  12  0 

1  -6  13  12  -30  0 

0  0  0  0  0  -1  . 


(2.36) 


(2.37) 
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The  next  task  is  to  show  that  M  =  |(M  +  M^)  is  negative  definite.  We  write  the  symmetric 
matrix  M  as  a,  sum  of  five  symmetric  matrices, 


M  —  [^0-^1  +  2M2  +  (24  —  ^o)M3  +  M4  +  Ms 


(2.38) 


We  shall  show  that  Mi  is  negative  definite,  and  that  Mj{j  =  2, ...  5)  are  non-positive  definite. 
The  M’s  are  given  by 


Ml  = 


1 

'2po 

0 

0 

0 

0 


0 

-2 

1 

0 

0 


0 

1 

-2 

1 

0 


0 

1 

-2 

1 


0 

0 

1 

-2 


1  -2  1 
1  -2 

0  0 


0 

0 

1 

■2/?o 


Mf  +  Mi  +  Mf  (2.39) 


where  M^  = 
middle  block. 


[-1/2^0  o' 
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[0  0  1 

0 

0 

CM 

1 

0 

_ 1 

and  Ml  is  the  remaining  {N—2)  x  {N—2) 


I 
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0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1 

1 

0 

0 

0 

0 

0 

1 

-2  1 

1  -2  1 

Mz  =  (2.41) 


1  -2  1 

1  -1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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-1/2 

0 

0 

0 

0 

0 

0 

-30  +  2^ 

2Q!2'’2 

12-/3 

-(o:2r3  +  0:372) 

13 

-(027-4  +  0:472) 

-6 

-(0275  +  0572) 

1 

0 

12 

-{0i2T^  +  azT2) 

-60  +  20 

—2otzTz 

32-0 

-(037-4  +  0473) 

-2 

-(0375  +  0:573) 

0 

0 

13 

-{a2T4  +  0472) 

32-0 

-(0:374  +  0:473) 

-60  +  20 

—20:474 

32-0 

-(0475  +  a574) 

-2 

0 

-6 

-{(X2T5  +  0:572) 

-2 

-(o:3r5  +  0:573) 

32-0 

-(a4T5  +  0:574) 

-58  +  ^ 

-2a57-5 

28-/3 

mi 

1 

0 

-2 

23-0 

-26  +  ^ 

(2.42) 
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Ms  = 


'  -26  +  0 

1 

oo 

-2 

0 

1 

0 

28-0 

-58  +  20 

—2aN-4TN-4 

82-0 
—  {ocN-3TN-4 
+aN-4TN-3) 

-2 

—{0CN-2TN-4 

+Ot]^-4TN-2) 

-6 

—  {cxN-1TN-4 
+ajV-4TAr-l) 

0 

-2 

82-0 
—  iotN~3TN-4 
+0CN-4:TN-3 

-60  -f  20 
—2aN-3TN-3 

82-0 

—{0CN-2TN-3 

+OCN-3TN-2) 

13 

—  {o‘N-l'rN-3 
+<^JV-3TZV-i) 

0 

0 

-2 

—  (0CN-2TN-4 
+aN-4TN-2) 

82-0 
—  {<y-N-2TN-3 
+aN-zTN-2) 

-60  -K  20 
—2aN-2TN-2 

12-0 
—  {oiN-lTN-2 
+o:N-2TN-i) 

0 

1 

-6 

—  {oiN-iTN-4 
d-Q;Ar-477v-i) 

13 

—  {o^N-1TN-3 
+OCN-3TN-1 

12-0 
—  {o^N-iT‘N-2 
+aN-2TN-l) 

! 

-30  +  20 
—2oin-iTn-i 

0 

Let  us  consider  Mi  -  see  (2.39);  it  may  be  decomposed  as  follows: 


(2.43) 


1  -1 


Mi  =  - 


(2.44) 


-1  1 


The  last  matrix  in  non-positive  definite.  The  first  term  is  a  product  of  a  regular  matrix  with 
its  transpose,  hence  its  negative  is  a  negative  definite  matrix.  Thus  we  established  that  Mi 
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is  negative  definite  for  any  finite  dimension  N.  All  its  eigenvalues  are  negative.  It  remains 
to  show  that  the  eigenvalues  of  (see  (2.38)  are  bounded  away  from  zero  by  a  constant 

as  h  0  {N  — >•  oo). 

Consider  a  symmetric  tridiagonal  matrix  S  with,  like  Mi,  constant  diagonals: 

b  a  0 
aba 

0  a  b  a  ak\ 

S=  .  .  (2.45) 


aba 
a  b 


Designate  by  Dj  the  determinant  of  the  upper-left  j  x  j  sub-matrix.  Thus  Di  =  6,  I>2  — 

det  ^  ?  ,  etc. 

a  b 

We  have  then  Di  =  6,  £>2  =  6^  -  and  in  general 


Dj  =  bDj-i  -  a?Dj-2 


(2.46) 


It  can  be  shown  (see  Appendix  I)  that  the  solution  to  the  recursion  relation  (2.46)  is 

Di  =  -4  [4  +  4]  P.47) 

«  L/^i  iA\ 

where 


/^i  =  —  [fe-fV62-4a2 


(2.48) 


(2.49) 
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A 


Ml  ~  M'i 


(2.50) 


B  =  — - [iD2-bDr)fi2  +  Dr]  (2.51) 

Ml  ~  M2 

We  have  already  shown  that  Mi  is  negative  definite.  The  eigenvalue  of  Mi  are  found  from 

det(M,  -  IX)  =  (-^  -  a)  •  det(iy',  -  XI)  ■  -  A^  =  0  (2.52) 

thus  either  A  =  — l/2/?o  <  0  (because  /?o  will  be  taken  positive)  or  A  =  eigenvalue  of  Mi  <  0. 
We  would  like  to  investigate  the  behavior  of  the  eigenvalues  of  -^pMi.  In  particular  we 
would  like  to  show  that  these  eigenvalues  (which  are  negative)  are  bounded  away  from  zero. 
To  show  this  we  analyze  the  behavior  of  Mi  —  A/  as  iV  increases.  We  now  take  S  =  Mi  —  XL 
Its  determinant  is  given  by  Dn-2-  Substituting  (2.48)-(2.51)  into  (2.47)  with  j  —  N  —  2  we 
get  after  some  elementary  manipulations 

=  (2.53) 

where 


p  =  ■v/4_ft2.  j)  —  _2  —  X-,  a  =  l 

r  =  yJlP  +  =  2 

6  —  tan“^(p/6) 


From  (2.52)  we  require 


Dn-2  —  0 


(2.54) 


(2.55) 
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This  is  equivalent,  see  (2.53),  to  requiring 


N-V 

From  the  definition  of  0  and  (2.54)  we  obtain 


k  =  l,...,N -2. 


N-1 


-A(A  +  4) 
2  -f-  A 


(A  <  0). 


Squaring  (2.57)  we  get  a  quadratic  equation  for  A,  the  solution  of  which  is 


A  =  —2  1  ±  1  +  tan 


N-1 


(2.56) 


(2.57) 


=  — 2  1  ±  cos 


N-1 


For  any  fixed  N,  the  smallest  values  of  |A|  is  given  by  (2.58)  for  jfc  =  1, 


-^max  —  ^  ^  \iV - 1 


As  N  increases,  we  have 


^  ^  2{N  -  1)2  (iV4 


'{N-iy 


-7r^h\ 


(2.58) 


(2.59) 


(2.60) 


Thus  the  eigenvalues  of  Mi/2ih^  (and  hence  of  Mij24:h?)  are  bounded  away  from  zero  by 
the  value  —  j . 

We  now  consider  M2-  One  can  verify  that 


M2  =  -M2MJ 


(2.61) 
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where 


0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0  0 

0  0  0  0  0  0  0 

0  0  0  0  0  0  1  0 

000000  -2  1  0 
0  1-21 
1  -2  1 

(2.62) 

0 

1  0 

2  0 

1  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

Therefore  M2  is  non-positive  definite.  In  a  similar  fashion  M3  is  non-positive  definite  because 

M3  =  -M3M3  (2.63) 


M2  = 

1 

0  1  -2 

0  0  1- 

0 


20 


with 


0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0  -1 

1  -1 

(2.64) 

1  -1 

1  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 

The  matrices  M4  and  are  N  x  N  matrices  with  zero  entries  except  for  6  x  6  upper-left 
(lower-right)  blocks.  It  is  sufficient  to  show  that  these  blocks  are  negative  definite.  This 
was  done  symbolically  using  the  Mathematica  software  and  plotted  for  0  <  7i,,7h  <  1  and 
/3o  =  1-  M4  and  M5  are  indeed  negative  definite  for,  0  <  <  1-  Thus  we  have  shown 

that  M  =  \{M  4-  M^)  is  indeed  negative  definite,  and  its  eigenvalues  are  bounded  away 
from  zero  by  (— 7r^/24),  even  as  N  00,  and  the  error  estimate  (2.13)  is  valid. 

3  The  Two  Dimensional  Case 


We  consider  the  inhomogeneous  diffusion  equation,  with  constant  coefficients,  in  a  domain 
n.  To  begin  with  we  shall  assume  that  0  is  convex  and  has  a  boundary  curve  dCl  G  C^. 
The  convexity  restriction  is  for  the  sake  of  simplicity  in  presenting  the  basic  idea;  it  will  be 
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removed  later.  We  thus  have 


dt  ' 

( &^u  d^u\ 

1  5x2  +  5^,2) f  >  0;  k>0 

(3.1a) 

u{x,y,0)  =  uo(x,y) 

(3.1b) 

u{x,y,t)\sci  = 

(3.1c) 

We  shall  refer  to  the  following  grid  representation: 


Figure  2:  Two  dimensional  grid. 


We  have  Mr  rows  and  Me  columns  inside  Cl.  Each  row  and  each  column  has  a  discreitized 
structure  as  in  the  one  1-D  case,  see  figure  1.  Let  the  number  of  grid  points  in  the  row 
be  denoted  by  Rk  and  similarly  let  the  number  of  grid  points  in  the  j*'^  column  be  Cj.  Let 
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the  solution  projection  be  designated  by  Uj,k{t).  By  U(t)  we  mean,  by  analogy  to  the  1-D 


U(t)  —  (ui,i,  U2,1,  .  .  . ,  Wl,2,  U2,2,  •  •  • ,  UR2,2'i  *  •  •  i  “I.Mb,  U2,Mr,  •  •  •  Ur^^^Mr) 

=  (ui,U2,...,umJ  (3.2) 

Thus,  we  have  arranged  the  solution  projection  array  in  vectors  according  to  rows,  starting 
from  the  bottom  of  ft. 

If  we  arrange  this  array  by  columns  (instead  of  rows)  we  will  have  the  following  structure 
u(‘')(t)  =  (ui,i,  Ui,2,  .  .  •  ,  W2,l5  “2,2,  •  •  • ,  “2,C2;  •  •  •  i  “Md,  “Me, 2,  •  •  •  ,  “McCji^J 

=  (uS"\u^"\...,uSj  (3.3) 

Since  U(‘=l(t)  is  just  a  permutation  of  U(t),  there  must  exist  an  orthogonal  matrix  P  such 


U(">(t)  =  PU  (3.4) 

If  the  length  of  U(t)  is  £,  then  P  is  an  ^  x  ^  matrix  whose  each  row  contains  i  —  1  zeros  and 
a  single  1  somewhere. 

The  second  derivative  operator  in  (3.1a)  is  represented  on  the  row  by  the 

differentiation  matrix  whose  structure  is  given  by  (2.14).  Similarly  let  d^fdy'^  be  given 
on  the  column  by  Df\  whose  structure  is  also  given  by  (2.14).  With  this  notation  the 
Laplacian  of  the  solution  projection  is: 

“b(^)  =  +  tW  + 
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(3.5) 


where 


)  1  r 

-v^y)  = 

n(x) 

J  L 

where  and  T>^y^  are  (£  x  £)  matrices  and  have  the  block  structures  shown.  and 
are  the  truncation  errors  associated  with  and  'D^y\  respectively.  We  now  call  attention 
to  the  fact  that  and  do  not  operate  on  the  same  vector.  This  is  fixed  using  (3.4): 

=  V^U  =  (PW  +  P^D(s')p)U  +  +  p'^T^y)  (3.7) 

Thus  (3.1a)  becomes,  by  analogy  to  (2.5), 

rIJ  T 

—  =  +  P^D(J')p)U  +  jfc(Tf)  +  P^ri^))  +  f(t)  (3.8) 

where  f (t)  is  f{x,  y;  t)  arranged  by  rows  as  a  vector. 

Before  proceeding  to  the  semi-discrete  problem  let  us  define: 

A/i*>  =  flW  -  (3.9) 

where  rx,,, ,  Al^  are  the  tl  and  Ai,  defined  in  section  2,  appropriate  to  the  row;  similarly 
for  Tji^  and  Ar^.  .  In  the  same  way,  define 

—  TBjAi,-  —  ttjAr.  (3.10) 

where  B  and  T  stand  for  bottom  and  top. 

We  can  now  write  the  semi-discrete  problem  by  analogy  to  (2.6) 

^  -h  P^A4(^)P)V  -I-  -f  +  f(t)  (3.11) 
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where  V  is  the  numerical  approximation  to  U; 


■  1  r 

;  (3.12) 

MiSi  J  [  MjSl  . 

and 

^  d"  TftiSila)?  •  •  •  j  ('TLfcS-t'*  d”  '^RkSRk}^  •  •  •  5  (TLmtSLmj^  d"  j 

=  [(rsigBi  d-  •  • . ,  {rSjgBj  +  tt^Etj),  •  •  • ,  {tbm.Sblmc  +  '^Tm.STmJ]  •  (3-13) 

Subtracting  (3.11)  from  (3.8)  we  get  in  a  fashion  similar  to  the  derivation  of  (2.9): 

dE 

—  =  +  kT2  (3.14) 

where  E  =  U  —  V  is  the  two  dimensional  array  of  the  errors,  €ij,  arranged  by  rows  as  a 
vector.  T2  is  proportional  to  the  truncation  error. 

The  time  change  of  ||  E  |p  is  given  by 

~  II  E  lp=  k{E,  (A4(")  +  P^A4(")P)E)  +  jfc(E,  T^)  (3.15) 

The  symmetric  part  of  +  P^M^^^P  is  given  by 

+  Af  W"')  +  M  (^)'’)P]  (3.16) 

Clearly  and  are  block-diagonal  matrices  with  typical  blocks 

given  by  and  .  We  have  already  shown  in  the  one  dimensional 

case  that  each  one  of  those  blocks  is  negative  definite  and  bounded  away  from  zero  by  7r^/24. 
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Therefore  the  operator  (3.16)  is  also  negative  definite  and  bounded  away  from  zero.  The 
rest  of  the  proof  follows  the  one  dimensional  case  and  thus  the  norm  of  the  error,  ||  £'  ||,  is 
bounded  by  a  constant. 

If  the  domain  Cl  is  not  convex  or  simply  connected  then  either  rows  or  columns,  or  both, 
may  be  “interrupted”  by  39,.  In  that  case  the  values  of  the  solution  on  each  “internal” 
interval  (see  figure  [3]  below)  are  taken  as  separate  vectors. 


Decomposing  “interrupted”  vectors  in  this  fashion  leaves  the  previous  analysis  unchanged. 
The  length  of  U  (or  is  again  £,  where  £  is  the  number  of  grid  nodes  inside  Cl.  The 
differentiation  and  permutation  matrices  remain  £x£.  Note  that  adding  more  “holes”  inside 
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dil  does  not  change  the  general  approach. 


4  Numerical  Example 

In  this  section  we  describe  numerical  results  for  the  following  problem: 

—  =  k{u^a>  +  Uyy)  +  f{x,  y, t),  {x,  y)en,  t  >  0,  (4.1) 

at 

where  0,  is  the  region  contained  between  a  circle  of  radius  tq  =  1/2  and  inner  circle  of  radius 
Vi  <  0.1.  The  inner  circle  is  not  concentric  with  the  outer  one.  Specifically  is  described  by 

{(a:  -  .5)"*  +  (y  -  .5)"*  <  1/4}  n  {(a:  -  .6)"  +  (y  -  .5)^  >  (.1  -  8f-  0  <  ^  <  .l)  (4.2) 

The  cartesian  grid  in  which  0  is  embedded  spans  0  <  a:,  y  <  1.  We  took  Ax  =  Ay,  and  ran 
several  cases  with  Ax  =  1/50,  1/75,  1/100.  The  geometry  thus  looks  as  follows: 

1 

OJ 

y=0 

x=0  OJ  0.6  1 

Figure  4: 
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The  source  function  f{x,y,t)  was  chosen  different  from  zero  so  that  we  could  assign  an 
exact  analytic  solution  to  (4.1).  This  enables  one  to  compute  the  error  Eij  =  Uij  —  Vij 
“exactly”  (to  machine  accuracy).  We  chose  k  =  1  and 

u{x,y,t)  =  1  -t-  cos(10f  —  lOx^  —  lOy^)  (4.3) 


This  leads  to 


f{x,y,t)  =  400(a:^  +  j/^)  cos(10f  —  lOx^  —  lOy^) 

—  50sin(10<  —  lOx^  —  lOy^)  (4.4) 

From  the  expression  for  u(x,y,t)  one  obtains  the  boundary  and  initial  conditions. 

The  problem  (4.1),  (4.2),  (4.4)  was  solved  using  both  a  “standard”  fourth  order  algorithm 
(a  2-D  version  of  (2.3))  and  the  new  “SAT,”  or  “bounded  error,”  approach  described  in 
Section  3.  The  temporal  advance  was  via  a  fourth  order  Runge-Kutta. 

The  standard  algorithm  was  run  for  Ax  =  1/50  and  a  range  of  0  <  5  <  .01  (.09  <  Vi  <  .1). 
We  found  that  for  S  >  .0017323,  the  runs  were  stable  and  the  error  bounded  for  “long”  times 
(10®  time  steps,  or  equivalently  t  =  2).  For  0  <  ^  <  .0017233  the  results  began  to  diverge 
exponentially  from  the  analytic  solution.  The  “point  of  departure”  depended  on  6.  A 
discussion  of  these  results  is  deferred  to  the  next  section.  Figures  5,6,7  show  the  T2-norm  of 
the  error  vs.  time  for  different  radii  of  the  inner  “hole.” 

The  same  configurations  were  also  run  using  the  “bounded  error”  algorithm  described  in 
Section  3  (see  eq.  (3.5)),  and  the  results  are  shown  in  figures  8,9,10,11.  It  is  seen  that  for 
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6^s  for  which  the  standard  methods  fails,  the  new  algorithm  still  has  a  bounded  error,  as 
predicted  by  the  theory. 

To  check  on  the  order  of  accuracy,  the  “SAT”  runs  (with  5  =  0)  were  repeated  for 
Ax  =  Ay  =  1/75  and  1/100.  Figure  12,13,  and  14  show  the  logarithmic  slope  of  the  2/2,  Ti 
and  Loo  errors  to  be  less  than  —4;  i.e.,  we  indeed  have  a  4‘^  order  method.  That  the  slopes 
are  larger  in  magnitude  than  4.5  is  attributed  to  the  fact  that  a.s  Arc  =  Ay  decreases  the 
percentage  of  “internal”  points  increases  (the  boundary  points  have  formally  only  3^“^  oder 
accuracy).  It  is  therefore  possible  that  if  the  number  of  grid  points  was  increased  much 
further,  the  slope  would  tend  to  —4.  Lack  of  computer  resources  prevented  checking  this 
point  further.  (For  Arc  =  0.01,  running  20,000  time  steps,  t  =  .1,  cpu  time  on  a  CRAY  YMP 
is  about  5  hours).  It  should  also  be  noted  that  the  “bounded-error”  algorithm  was  run  with 
a  time  step.  At,  twice  as  large  as  the  one  used  in  the  standard  scheme.  At  this  larger  At 
the  standard  scheme  “explodes”  immediately. 


err 


err 
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0.00025 
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0.00015 

0.0001 

0.00005 
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t 


Figure  5:  S  ^  0.0017325,  Standard 
scheme 


Figure  6:  ^  =  0.0017323,  Standard 
scheme 
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Figure  1:  8  =  0.0015,  Standard  scheme 
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Figure  8:  <5  =  0,  SAT  scheme 
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Figure  9:  5  =  0.0015,  SAT  scheme 


Figure  10:  8  =  0.0017323,  SAT  scheme 


Figure  11:  <5  =  0.0017325,  SAT  scheme 


Figure  12:  Order  of  accuracy  L\ 
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Figure  13:  Order  of  accuracy  Figure  14:  Order  of  accuracy  Loo 

A  study  of  the  effect  of  size  of  At  shows  that  the  instabilities  exhibited  above  are  due  to 
the  time-step  being  near  the  C.F.L.-limit.  It  is  interesting  that  this  C.F.L.-limit  depends  so 
strongly  on  the  geometry. 

5  Conclusions 

(i)  The  theoretical  results  show  that  one  has  to  be  very  careful  when  using  an  algorithm 
whose  differentiation  matrix,  or  rather  its  symmetric  part,  is  not  negative  definite.  For 
some  problems,  such  “standard”  schemes  will  give  good  answers  (i.e.,  bounded  errors) 
and  for  others  instability  will  set  in.  Thus,  for  example,  the  “standard”  scheme  for 
the  1-D  case  has  a  matrix  which,  for  all  0  <  7L,7ij  <  1,  though  not  negative  definite 
has  eigenvalues  with  negative  real  parts.  This  assures,  in  the  1-D  case,  the  temporally 
asymptotic  stability.  In  the  2-D  case,  even  though  each  of  the  block  sub-matrices  of 
the  i  X  i  x-and-t/  differentiation  matrices  has  only  negative  (real-part)  eigenvalues,  it 
is  not  assured  that  the  sum  of  the  two  £  x  £  matrices  will  have  this  property.  This 
depends,  among  other  things,  on  the  shape  of  the  domain  and  the  mesh  size  (because 


31 


the  mesh  size  determines,  for  a  given  geometry,  the  and  'Jrs  along  the  boundaries). 

Thus  that  we  might  have  the  “paradoxical”  situation,  that  for  a  given  domain  shape, 
successive  mesh  refinement  could  lead  to  instability  due  to  the  occurrence  of  destabi¬ 
lizing  7’s.  This  cannot  happen  if  one  constructs,  as  was  done  here,  a  scheme  whose 
differentiation  matrices  have  symmetric  parts  that  are  negative  definite. 

It  is  also  interesting  to  note  that  if  one  uses  explicit  standard  method  then  the  allow¬ 
able  C.F.L.  may  decrease  extremely  rapidly  with  change  in  the  geometry  that  causes 
decrease  in  the  7’s.  This  point  is  brought  out  in  figures  5  to  7. 

(ii)  Note  that  the  construction  of  the  2-D  algorithm,  and  its  analysis,  which  were  based 
on  the  1-D  case,  can  be  extended  in  a  similar  (albeit  more  complex)  fashion  to  higher 
dimensions. 

(iii)  Also  note  that  if  the  diffusion  coefficient  k,  in  the  equation 

Ut  =  kA^u 

is  a  function  of  the  spatial  coordinates,  k  =  k{x,y,z),  the  previous  analysis  goes 
through  but  the  energy  estimate  for  the  error  is  now  for  a  different,  but  equivalent 
norm. 
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Appendix  I 


We  start  with 

Dj  =  bDj-i  —  a^Dj-2 

with 

Di=b,D2  =  9-0^ 

We  associate  with  (A.l)  a  generating  function  /(x), 

C» 

fi^)  =  Dj+tx^ 
j=o 

Multiplying  (A.l)  by  for  each  j  >  3,  and  summing  both  sides  we  obtain: 

f-Di-  Da,x  -Dx  _2  X 

2  —  ^  ^  J 

nr* 


leading  to 


1  Dx  d"  {D2  —  bDx'jx 
—  (6/o^)x  +  (l/a^)_ 


(A.1) 

(A.2) 

(A.3) 

(A.4) 


1  Dx  +  {D2  —  bDx)x 
a?  [x  —  Ui)(a:  —  U2) 


where  ui,U2  are  given  by  (2.48),  (2.49). 
We  may  also  present  /  by 


—  ^  +  jg 

a?  [x  —  Ux)  (x  —  U2) 


(A.5) 


(A.6) 
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Comparing  (A.6)  and  (A.5)  we  get  expression  for  A  and  B  as  given  in  (2.50),  (2.51).  Ex¬ 
panding  the  denominator  in  (A.6)  we  get  the  following  series  for  / 


,  1  ^  /  A  B 

"  ■>=  §  ■ 


(A.7) 


from  which  it  immediately  follows  (see  (A. 3))  that 


r.  1  [A  B 

—  — 2  I  ~7  — 7 
\ul  U2 


(A.8) 
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